#######################################################################################################################
############ Governmental Responses to Terrorism in Autocracies: Evidence from China ##################################
############ Online Appendix ##########################################################################################
## Philip B.K. Potter, Associate Professor, Department of Politics, University of Virginia
## Chen Wang, PhD Candidate, Department of Politics, University of Virginia
#######################################################################################################################

## Load packages ######################################################################################################
library(expss)
library(foreign)
library(survival)
library(KMsurv)
library(MASS)
library(stargazer)
library(pastecs)
library(dplyr)
library(tidyr)
library(rmarkdown)
library(knitr)
library(ggplot2)
library(survminer)
library(pastecs)
library(pscl)
library(simPH)
library(ggthemes)
library(margins)
library(lubridate)
library(sandwich)
library(lmtest)
library(gridExtra)
######################################################################################################################
## Load data #########################################################################################################

## Set working directory
setwd()

## Load data
d1<- readRDS("data_governmental_response_terrorism.RDS")

## Check variable names and labels
names(d1)

checkvariable <- data.frame(matrix(nrow = 32, ncol = 2))
colnames(checkvariable) <- c("variable name", "labels")
for (x in 1:32) {
  checkvariable[x,1]<- colnames(d1[x])
  checkvariable[x,2]<- var_lab(d1[,x])
}
View(checkvariable)

######################################################################################################################

## Generate Appendix Figure A.1 ######################################################################################
names(d1)
d1$other<- d1$incendiary+d1$gun+d1$vehicle+d1$poison

dweap<- select(d1, year, knive, bombing, other, weap_unknown)
dweap<- group_by(dweap, year)
dweap<- summarise(dweap, knive=sum(knive), bombing=sum(bombing),
                  other=sum(other), unknown=sum(weap_unknown))
dweap_ts<- data.frame("year"=1990:2014)
dweap_ts<- apply_labels(dweap_ts, year="Year of the event")
dweap_ts<- left_join(dweap_ts, dweap, by="year")
dweap_ts$knive<- ifelse(is.na(dweap_ts$knive)==TRUE,0,dweap_ts$knive)
dweap_ts$bombing<- ifelse(is.na(dweap_ts$bombing)==TRUE,0,dweap_ts$bombing)
dweap_ts$other<- ifelse(is.na(dweap_ts$other)==TRUE,0,dweap_ts$other)
dweap_ts$unknown<- ifelse(is.na(dweap_ts$unknown)==TRUE,0,dweap_ts$unknown)

bplot<- matrix(rbind(as.vector(t(dweap_ts[2])),as.vector(t(dweap_ts[3])),
                  as.vector(t(dweap_ts[4])),as.vector(t(dweap_ts[5]))),
                 nrow=4,
                 ncol=25)

par(mar=c(5,6,4,2)+0.2)
barplot(bplot,beside = FALSE, density = c(30,30,30,30),angle = c(45,0,120,90),xlab = "Year",ylab = "Count",
        legend.text = c("Knives","Bombing","Unknown","Other"),args.legend = list(x = 20,y=30,cex = 1.2),
        names.arg = c(1990:2014),cex.axis =2,cex.name=2,cex.lab = 2)
######################################################################################################################

## Generate Appendix Figure A.2 ######################################################################################
dreport<- subset(d1,delta==1)
par(mar=c(5,6,4,2)+0.2,mfrow=c(1,2))
hist(dreport$duration, xlab="Number of Days", breaks = 20,
     main="Histogram of Reporting Duration" , freq=TRUE, cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2)
sfit <- survfit(Surv(duration, delta)~1, data=d1)
plot(sfit, xlim=c(0, 30),main="Non-Parametric Model: Kaplan-Meier Estimator ",
     ylab = "Probability of Non-reporting",xlab = "Number of Days",
     cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2)
######################################################################################################################

## PH Tests: Appendix Table A.4 ######################################################################################
d1$cgrowth<- d1$growth-mean(d1$growth)                   # centered GDP growth
d1$cmcci<- d1$cci-mean(d1$cci)                           # centered Consumer Confidence Index
d1$cliindex<- d1$liindex-mean(d1$liindex)                # centered Li-index
d1$cmajor<- d1$majority_freq-mean(d1$majority_freq)      # centered Majority Frequency
d1$ctwe<- d1$majority_g20-mean(d1$majority_g20)          # centered Majority Frequency among G20 countries
d1$csc<- d1$majority_scmember-mean(d1$majority_scmember) # centered Majority Frequency among Security Council Members


cox1<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+cluster(id), data = d1,method="breslow")
cox2<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
               internet_ratio+sensitive+cluster(id), data = d1,method="breslow")
cox3<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
cox4<- coxph(Surv(duration, delta)~ cliindex+cmajor+cliindex:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
cox5<- coxph(Surv(duration, delta)~ cmcci+cmajor+cmcci:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
cox6<- coxph(Surv(duration, delta)~ cgrowth+ctwe+cgrowth:ctwe+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
cox7<- coxph(Surv(duration, delta)~ cgrowth+csc+cgrowth:csc+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")

d1$cnatural<- (d1$natural_disaster_days)-mean((d1$natural_disaster_days))       # centered Natural Disaster
d1$cuschina<- d1$uschinadis-mean(d1$uschinadis)                                 # centered US-China Distance
d1$cexternal<- d1$external_condition-mean(d1$external_condition, na.rm = TRUE)  # centered external condition
d1$cdomestic<- d1$domestic_condition-mean(d1$domestic_condition, na.rm=TRUE)    # centered domestic condition

cox8<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+idp+success+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
cox9<- coxph(Surv(duration, delta)~ cgrowth+cuschina+cgrowth:cuschina+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
cox10<- coxph(Surv(duration, delta)~ cnatural+cmajor+cnatural:cmajor+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
cox11<- coxph(Surv(duration, delta)~ cnatural+cuschina+cnatural:cuschina+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
cox12<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1[which(d1$year>1994),],method="breslow")
cox13<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+growth+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1[which(d1$year>2000),],method="breslow")

ph1<-cox.zph(cox1, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph2<-cox.zph(cox2, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph3<-cox.zph(cox3, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph4<-cox.zph(cox4, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph5<-cox.zph(cox5, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph6<-cox.zph(cox6, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph7<-cox.zph(cox7, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph8<-cox.zph(cox8, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph9<-cox.zph(cox9, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph10<-cox.zph(cox10, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph11<-cox.zph(cox11, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph12<-cox.zph(cox12, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)
ph13<-cox.zph(cox13, transform="identity", terms=TRUE, singledf=FALSE, global=TRUE)

## Get the Global Test for each model
phtest<- data.frame("Model"=c("Model 1", "Model 2","Model 3","Model 4","Model 5","Model 6",
                              "Model 7","Model 8","Model 9","Model 10","Model 11","Model 12",
                              "Model 13"),
                    rbind(ph1$table[4,], ph2$table[6,], ph3$table[10,],ph4$table[10,],
                          ph5$table[10,],ph6$table[10,],ph7$table[10,],ph8$table[12,],
                          ph9$table[10,],ph10$table[11,],ph11$table[11,],ph12$table[11,],
                          ph13$table[11,])
)
library(xtable)
xtable(phtest)
######################################################################################################################

## Generate Figure A.3 ###############################################################################################

## Get the plot for Model 6 in Table 1
cox6<- coxph(Surv(duration, delta)~ cgrowth+ctwe+cgrowth:ctwe+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id),
             data = d1,method="breslow")
summary(cox6)

pdata <- data.frame(cgrowth=c(mean(d1$cgrowth,na.rm=TRUE)+0.5*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-0.5*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)+0.5*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-0.5*sd(d1$cgrowth,na.rm=TRUE)),
                     ctwe=c(mean(d1$ctwe,na.rm=TRUE)-sd(d1$ctwe,na.rm=TRUE),
                            mean(d1$ctwe,na.rm=TRUE)-sd(d1$ctwe,na.rm=TRUE),
                            mean(d1$ctwe,na.rm=TRUE)+sd(d1$ctwe,na.rm=TRUE),
                            mean(d1$ctwe,na.rm=TRUE)+sd(d1$ctwe,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing)
                     
) 

fit0<- survival::survfit(cox6, newdata=pdata,se.fit=TRUE, conf.int=0.95,
                         conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("High Growth & Majority: Low Frequency",
                                 "Low Growth & Majority: Low Frequency",
                                 "High Growth & Majority: High Frequency",
                                 "Low Growth & Majority: High Frequency"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))

pcox6<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                          labels=c("High Growth \nHigh Majority Frequency",
                                                   "High Growth \nLow Majority Frequency",
                                                   "Low Growth \nHigh Majority Frequency",
                                                   "Low Growth \nLow Majority Frequency"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 6: G20 Majority Frequency")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = -6)))+
  labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=26),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.key.height=unit(4,"line"),
        legend.position="right",legend.justification = c("left"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(hjust = 0, vjust=-0.6,face = "italic"))


pcox6

## Get the plot for Model 7 in Table 1
cox7<- coxph(Surv(duration, delta)~ cgrowth+csc+cgrowth:csc+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")

pdata <- data.frame(cgrowth=c(mean(d1$cgrowth,na.rm=TRUE)+1.5*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-1*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)+1.5*sd(d1$cgrowth,na.rm=TRUE),
                               mean(d1$cgrowth,na.rm=TRUE)-1*sd(d1$cgrowth,na.rm=TRUE)),
                     csc=c(mean(d1$csc,na.rm=TRUE)-1*sd(d1$csc,na.rm=TRUE),
                           mean(d1$csc,na.rm=TRUE)-1*sd(d1$csc,na.rm=TRUE),
                           mean(d1$csc,na.rm=TRUE)+1.5*sd(d1$csc,na.rm=TRUE),
                           mean(d1$csc,na.rm=TRUE)+1.5*sd(d1$csc,na.rm=TRUE)),
                     internet_ratio=mean(d1$internet_ratio),
                     casualty=mean(d1$casualty,na.rm=TRUE),
                     sensitive=mean(d1$sensitive),
                     urban=mean(d1$urban),
                     target_civilian=mean(d1$target_civilian),
                     bombing=mean(d1$bombing)
                     
) 


fit0<- survival::survfit(cox7, newdata=pdata,se.fit=TRUE, conf.int=0.95,
                         conf.type=c("plain"),censor=FALSE, type="kaplan-meier")

survdata<- data.frame(Days=rep(c(1,2,3,4,19,21), 4),
                      Case=rep(c("High Growth & Majority: Low Frequency",
                                 "Low Growth & Majority: Low Frequency",
                                 "High Growth & Majority: High Frequency",
                                 "Low Growth & Majority: High Frequency"), each=6),
                      surv=c(fit0$surv[,1:4]),
                      upper=c(fit0$upper[,1:4]),
                      lower=c(fit0$lower[,1:4]))



pcox7<- ggplot(survdata,aes(x=Days,y=surv,linetype=Case))+
  geom_line(size=1)+scale_linetype_manual(values=c("solid","dashed","dotted","twodash"),
                                          labels=c("High Growth \nHigh Majority Frequency",
                                                   "High Growth \nLow Majority Frequency",
                                                   "Low Growth \nHigh Majority Frequency",
                                                   "Low Growth \nLow Majority Frequency"))+
  xlab("Number of days after the attack")+
  ylab("Probability of non-reporting")+ggtitle("Model 7: Security Council Majority Frequency")+
  #geom_hline(yintercept=0, linetype="dotted", color = "red",size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2, fill="grey30")+
  theme(axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = -6)))+
  #labs(caption = "Note: Shaded area is the 95% confidence interval.", face = "italic")+
  theme(panel.background = element_rect(fill = "white", colour = "black"), 
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        text = element_text(size=26),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.8, linetype="solid", colour ="black"),
        legend.text=element_text(size=18),
        legend.key.height=unit(4,"line"),
        legend.position="right",legend.justification = c("left"),
        legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"))


pcox7

grid.arrange(pcox6,pcox7, nrow = 1)
######################################################################################################################

## Models in Appendix Table A.5 ######################################################################################

## Alternative Models without GDP Growth as an control
coxA1<- coxph(Surv(duration, delta)~ cnatural+cmajor+cnatural:cmajor+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
summary(coxA1)

coxA2<- coxph(Surv(duration, delta)~ cnatural+cuschina+cnatural:cuschina+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1,method="breslow")
summary(coxA2)

coxA3<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1[which(d1$year>1994),],method="breslow")
summary(coxA3)

coxA4<- coxph(Surv(duration, delta)~ cdomestic+cexternal+
                cdomestic:cexternal+
                internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
              data = d1[which(d1$year>2000),],method="breslow")
summary(coxA4)
######################################################################################################################

## Logit Models in Table A.6 #########################################################################################
lg1<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor,
          data=d1, family="binomial")
summary(lg1)
coeftest(lg1, vcov = vcovCL(lg1, cluster =d1$id)) 

lg2<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor+internet_ratio+sensitive,
          data=d1, family="binomial")
summary(lg2)
coeftest(lg2, vcov = vcovCL(lg2, cluster =d1$id)) 

lg3<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor+casualty+urban+target_civilian+bombing,
          data=d1, family="binomial")
summary(lg3)
coeftest(lg3, vcov = vcovCL(lg3, cluster =d1$id))

lg4<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor+internet_ratio+sensitive+casualty+urban+target_civilian+bombing,
          data=d1, family="binomial")
summary(lg4)
coeftest(lg4, vcov = vcovCL(lg4, cluster =d1$id)) 
######################################################################################################################

## Generate Figure A.4 ###############################################################################################
lg4<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor+internet_ratio+sensitive+casualty+urban+target_civilian+bombing,
          data=d1, family="binomial")
summary(lg4)
coeftest(lg4, vcov = vcovCL(lg4, cluster =d1$id)) 


# For Majority Low
cgrowthvalue <- seq(-3, 4.15, by=1)
logit.prob <- sapply(cgrowthvalue, FUN=function(x){
  mean(predict(lg4, type = "response", 
               newdata = mutate(d1, cgrowth = x,cmajor=mean(d1$cmajor)-1.5*sd(d1$cmajor))), na.rm=TRUE)
})

QI <- data.frame(cgrowth = cgrowthvalue,
                 cmajor=rep("Low Majority Frequency", 8),
                 logit.prob = logit.prob)

B <- coef(lg4)
V <- vcovCL(lg4, cluster =d1$id)

sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)


logit2 <- lg4
sim.qi <- apply(sim.coefs, 1, FUN=function(x){
  logit2$coefficients <- x
  logit.prob <- sapply(cgrowthvalue, FUN=function(y){
    mean(predict(logit2, type = "response", 
                 newdata = mutate(d1, cgrowth = y,cmajor=mean(d1$cmajor)-1.5*sd(d1$cmajor))), na.rm=TRUE)
  })
})
sim.qi <- t(sim.qi)

logit.prob.se <- apply(sim.qi, 2, sd)
logit.prob.lb <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .025)
})
logit.prob.ub <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .975)
})
QI <- mutate(QI, 
             logit.prob.se = logit.prob.se,
             logit.prob.lb = logit.prob.lb,
             logit.prob.ub = logit.prob.ub)


#For Majority High
logit.prob2 <- sapply(cgrowthvalue, FUN=function(x){
  mean(predict(lg4, type = "response", 
               newdata = mutate(d1, cgrowth = x,cmajor=mean(d1$cmajor)+1.5*sd(d1$cmajor))), na.rm=TRUE)
})

QI2 <- data.frame(cgrowth = cgrowthvalue,
                  cmajor=rep("High Majority Frequency", 8),
                  logit.prob = logit.prob2)

logit2 <- lg4
sim.qi <- apply(sim.coefs, 1, FUN=function(x){
  logit2$coefficients <- x
  logit.prob <- sapply(cgrowthvalue, FUN=function(y){
    mean(predict(logit2, type = "response", 
                 newdata = mutate(d1, cgrowth = y,cmajor=mean(d1$cmajor)+1.5*sd(d1$cmajor))), na.rm=TRUE)
  })
})
sim.qi <- t(sim.qi)

logit.prob.se <- apply(sim.qi, 2, sd)
logit.prob.lb <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .025)
})
logit.prob.ub <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .975)
})
QI2 <- mutate(QI2, 
              logit.prob.se = logit.prob.se,
              logit.prob.lb = logit.prob.lb,
              logit.prob.ub = logit.prob.ub)

QIcomb<- rbind(QI, QI2)

QIcomb$growth<- QIcomb$cgrowth+mean(d1$growth)

pA4<- ggplot(QIcomb,aes(x=growth,y=logit.prob,linetype=cmajor))+
  geom_line(size=1.2)+scale_linetype_manual(values=c("solid","dotted"))+
  xlab("GDP growth rate (non-centered value)")+
  ylab("Probability of being reported")+
  geom_ribbon(aes(ymin=logit.prob.lb,ymax=logit.prob.ub),alpha=0.2)+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 20, r = 30, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)),
        text = element_text(size=30),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=26),
        legend.position="right",legend.key.size=unit(1.2, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 18, vjust=-6))
pA4
annotate_figure(pA4,  
                bottom = text_grob("Note: Shaded areas are 95% confidence intervals generated from 1000 draws of new coefficients from the posterior.", 
                                   hjust = -0.1, vjust=0.2,x = 0, size = 14))
######################################################################################################################

## Generate Figure A.5 ###############################################################################################
## Hainmueller el al. (2019) diagonostic
library(interflex)

pA5<- inter.raw(Y = "delta", D = "growth", X = "majority_freq", data = d1, nbins = 2,
                    Ylabel = "Probability of being reported", Dlabel = "GDP Growth Rate", Xlabel="Majority Frequency", 
                    theme.bw = TRUE, show.grid = FALSE,cex.lab=1.4, cex.axis = 1.2)
pA5+theme(text = element_text(size=26))
######################################################################################################################

## Generate Figure A.6 ###############################################################################################
## Binning Estimator
g1<- inter.binning(data=d1, Y = "delta", D = "cgrowth", X = "cmajor", treat.type = "continuous", base = NULL,
                   Z = c("internet_ratio","sensitive","casualty","target_civilian","bombing","urban"), 
                   FE = NULL, weights = NULL, full.moderate = TRUE,
                   na.rm = TRUE, Xunif = FALSE, nbins = 3, cutoffs = NULL, CI = TRUE,
                   vartype = "cluster", nboots= 200, parallel = TRUE, cores =4,
                   cl = "id", time = NULL, pairwise = TRUE, wald = TRUE, predict = FALSE,
                   D.ref = NULL, figure = TRUE, order = NULL, subtitles = "Binning Estimator", 
                   show.subtitles = NULL, Xdistr = "histogram", main = NULL,
                   Ylabel = NULL, Dlabel = NULL, Xlabel = "continuous", xlab = "Majority Frequency (non-centered value)",
                   ylab = "Marginal effect of GDP growth rate\n on the probability of being reported", 
                   xlim = NULL, ylim = NULL, theme.bw = FALSE,
                   show.grid = TRUE, cex.main = NULL, cex.sub = NULL, cex.lab = 1.4, cex.axis = 1.2,
                   bin.labs = TRUE, interval = NULL, file = NULL, ncols = NULL, pool = FALSE,
                   color = NULL, jitter = FALSE, legend.title = NULL)
g1$graph+scale_x_continuous(limits = c(-23, 25.5, by=10),
                            breaks = c(-20, -10, 0, 10, 20),
                            labels = c(30,40,50,60,70))
gg1<-g1$graph+ggtitle("Binning Estimator")+theme(plot.title = element_text(size = 18))+
  scale_x_continuous(limits = c(-23, 25.5, by=10),
                     breaks = c(-20, -10, 0, 10, 20),
                     labels = c(30,40,50,60,70))
gg1

## Kernel Estimator
## Since the confidence intervals are bootstrapped, it may look slightly different after every draw.
g2<- inter.kernel(data=d1, Y = "delta", D = "cgrowth", X = "cmajor", treat.type = "continuous", base = NULL,
                  Z = c("internet_ratio","sensitive","casualty","target_civilian","bombing","urban"), 
                  FE = NULL, weights = NULL, full.moderate = TRUE,
                  na.rm = TRUE, Xunif = FALSE, CI = TRUE, conf.level = 0.95, cl = NULL,
                  CV.method =NULL, kfold = 10, grid = 30, neval = 50, nboots = 200, parallel = TRUE,
                  cores = 4, seed = 02139, bw = NULL, bw.adaptive = TRUE, quantile.eval = FALSE,
                  metric = "MSPE", predict = FALSE, D.ref = NULL, figure = TRUE, order = NULL,
                  subtitles = NULL, show.subtitles = NULL, Xdistr = "histogram", main = NULL,
                  Ylabel = NULL, Dlabel = NULL, Xlabel = NULL, xlab = "Majority Frequency (non-centered value)", 
                  ylab = "Marginal effect of GDP growth rate\n on the probability of being reported",
                  xlim = NULL, ylim = NULL, theme.bw = FALSE, interval = NULL, show.grid = TRUE,
                  cex.main = NULL, cex.sub = NULL, cex.lab = 1.4, cex.axis = 1.2,
                  file = NULL, ncols = NULL, pool = FALSE, color = NULL, legend.title = NULL,
                  diff.values = NULL, percentile = FALSE)

g2$graph+scale_x_continuous(limits = c(-23, 25.5, by=10),
                            breaks = c(-20, -10, 0, 10, 20),
                            labels = c(30,40,50,60,70))

gg2<- g2$graph+ggtitle("Kernel Estimator")+theme(plot.title = element_text(size = 18))+
  scale_x_continuous(limits = c(-23, 25.5, by=10),
                     breaks = c(-20, -10, 0, 10, 20),
                     labels = c(30,40,50,60,70))
gg2

## Logit Model
lg4<- glm(delta ~ cgrowth+cmajor+cgrowth:cmajor+internet_ratio+sensitive+casualty+urban+target_civilian+bombing,
          data=d1, family="binomial")
summary(lg4)
coeftest(lg4, vcov = vcovCL(lg4, cluster =d1$id)) 


y.star <- predict(lg4, type="link", newdata=mutate(d1, cmajor=10))
coef.growth <- coef(lg4)[["cgrowth"]]
coef.interaction <- coef(lg4)[["cgrowth:cmajor"]]
marg.prob.growth <- (coef.growth+coef.interaction*10)*exp(-y.star) / (1 + exp(-y.star))^2
mean(marg.prob.growth, na.rm=TRUE)

tenvalues <- seq(-20, 26, by=2)
logit2 <- lg4
B <- coef(lg4)
V <- vcovCL(lg4, cluster=d1$id)
sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)
growth.df <- sapply(tenvalues, FUN=function(majorityval){
  y.star <- predict(lg4, type="link", newdata=mutate(d1, cmajor=majorityval))
  coef.growth <- coef(lg4)[["cgrowth"]]
  coef.interaction <- coef(lg4)[["cgrowth:cmajor"]]
  marg.prob.growth <- (coef.growth+coef.interaction*majorityval)*exp(-y.star) / (1 + exp(-y.star))^2
  point.est <- mean(marg.prob.growth, na.rm=TRUE)
  marg.prob.growth.sim <- apply(sim.coefs, 1, FUN=function(x){
    logit2$coefficients <- x
    y.star <- predict(logit2, type="link", newdata=mutate(d1, cmajor=majorityval))
    coef.growth <- coef(lg4)[["cgrowth"]]
    coef.interaction <- coef(lg4)[["cgrowth:cmajor"]]
    marg.prob.growth <- (coef.growth+coef.interaction*majorityval)*exp(-y.star) / (1 + exp(-y.star))^2
    return(mean(marg.prob.growth, na.rm=TRUE))
  })
  se <- sd(marg.prob.growth.sim)
  lb <- quantile(marg.prob.growth.sim, .025)
  ub <- quantile(marg.prob.growth.sim, .975)
  return(c(point.est, se, lb, ub))
})
growth.df <- data.frame(tenvalues, t(growth.df))
colnames(growth.df) <- c("majorityval", "point.est", "se", "lb", "ub")

growth.df$major<- growth.df$majorityval+mean(d1$majority_freq)
g3<- ggplot(growth.df,aes(x=major,y=point.est))+
  geom_line(size=1.2)+scale_linetype_manual(values=c("solid"))+
  geom_hline(yintercept = 0, color="red", linetype="dotted")+
  ggtitle("Logit Model (Model A8)")+
  xlab("Majority Frequency (non-centered value)")+
  ylab("Marginal effect of GDP growth rate\n on the probability of being reported")+
  geom_ribbon(aes(ymin=lb,ymax=ub),alpha=0.2)+
  theme(plot.title = element_text(size=18),
        axis.title.y = element_text(margin = margin(t = 20, r = 30, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)),
        text = element_text(size=22),legend.title = element_blank(),
        legend.background = element_rect(fill="gray",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=26),
        legend.position="right",legend.key.size=unit(1.2, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 16, vjust=0, hjust=0))
g3

library(cowplot)
gg<- plot_grid(g3, gg1,gg2, nrow =1) # aligning vertically along the left axis
gg
annotate_figure(gg,  
                bottom = text_grob("Note: shaded areas present 95% confidence intervals", 
                                   hjust = -0.1, vjust=0.2,x = 0, size = 14))
######################################################################################################################

## Robustness check with EVC data ####################################################################################
drc<- readRDS("evc_robustness.RDS")
names(drc)


drc$cgrowth<- drc$growth-mean(drc$growth)                   # centered GDP growth
drc$cmajor<- drc$majority_freq-mean(drc$majority_freq)      # centered Majority Frequency


## The same as Model-1 in Table 1
coxA9<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+cluster(id), data = d1,method="breslow")
summary(coxA9)

coxA10<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+cluster(id), data = drc,method="breslow")
summary(coxA10)

## The same as Model-3 in Table 1
coxA11<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
               internet_ratio+sensitive+casualty+urban+target_civilian+bombing+cluster(id), 
             data = d1,method="breslow")
summary(coxA11)

coxA12<- coxph(Surv(duration, delta)~ cgrowth+cmajor+cgrowth:cmajor+
                 internet_ratio+sensitive+urban+target_civilian+bombing+cluster(id), 
               data = drc,method="breslow")
summary(coxA12)
